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ABSTRACT 


Numerical  tables  available  for  M/E, /c  queueing  systems 
are  discussed.   A  new  approximation  method  for  steady— state 
information  and  waiting  time  distribution  of  this  queueing 
system  was  developed.   Validity  of  approximation  was  estab- 
lished directly  for  the  large  waiting  times  and  by  simula- 
tion for  the  smaller  values.   The  developed  method  enables 
one  to  find  delay  probability,  expected  number  in  the 
queue  and  in  the  system,  expected  time  to  be  spent  in  the 
queue  and  in  the  system,  and  probability  of  waiting  for 
more  than  a  specified  time  t. 
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I.   INTRODUCTION  AND  SUMMARY 

Multi— server  queues  constitute  a  large  proportion  of 
queueing  systems  which  arise  in  practice.   Among  those, 
queues  with  Poisson  arrivals  and  Erlang  service  times  (which 
will  be  denoted  as  M/E,/c)  occupy  an  important  position,  as 
seen  in  the  banks,  airport  check— in  counters,  hotels,  super- 
markets etc.   However,  no  significant  theoretical  work  was 
done  until  the  late  60 's.   In  recent  years,  methods  for  the 
steady  state  solution  of  M/E, /c  queues  became  available,  fol- 
lowed by  the  numerical  tables  which  were  obtained  by  imple- 
mentation of  these  methods.   Even  so,  there  is  still  some 
computational  difficulty  involved  obtaining  some  particular 
information  about  above— mentioned  queues,  such  as  probabili- 
ty of  delay  exceeding  a  specified  time  length. 

The  importance  of  the  M/Ek/c  queue  is  due  to  the  fact 
that  it  is  used  to  model  any  queueing  system  whose  service 
time  distribution  is  believed  to  be  unimodal.   The  solutions 
are  known  for  extreme  values  of  k:  exponential  service  time 
distribution  for  k=l  and  constant  service  times  as  k-*-00. 
Usable  solutions  for  multi— server  queues  having  either  expo- 
nential or  constant  service  times  are  readily  available. 
The  M/E, /c  queue  is  also  important  by  providing  a  large 
variety  of  service  time  distributions  in  between  these  two 
extremes. 


The  purpose  of  this  study  is,  by  means  of  computer  simu- 
lation and  examination  of  numerical  tables  published  earlier, 
to  present  a  simple  and  accurate  computation  method  for  ob- 
taining information  about  steady— state  solution  of  the  M/E,  /c 
queue. 

For  this  purpose,  a  general  description  of  M/Ek/c  queue— 
ing  model  with  definitions  and  assumptions  for  the  analysis 
of  it  are  given  in  Chapter  II,  followed  by  a  summary  of  for- 
mer studies. 

In  Chapter  III,  the  numerical  tables  are  analyzed  and 
some  conclusions  are  drawn  in  terms  of  a  simple  form  approxi- 
mation for  steady— state  distribution  of  waiting  times. 

This  approximation  for  the  distribution  of  waiting  times 
needed  verification  for  small  values  of  them.   These  cases 
were  studied  by  computer  simulation.   The  procedure  and  the 
results  of  simulation  are  given  in  Chapter  IV. 

The  results  of  the  analysis  of  numerical  tables  and  simu- 
lation are  combined  and  generalized  in  Chapter  V.   Then  the 
computational  method  developed  by  this  generalization  is  de- 
scribed on  a  step— by— step  basis.   With  this  method,  there  is 
no  need  for  numerical  tables,  computations  are  very  simple, 
and  the  results  are  accurate  for  most  practical  purposes 
(the  error  being  about  2%  in  general) .   To  demonstrate  the 
advantages  of  the  proposed  method,  a  comparison  of  it  with 
the  other  method  which  uses  numerical  tables  is  given  at  the 
end  of  that  chapter. 
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II.   The  M/E,  /c  Queue 

The  description  of  the  M/E,/c  queueing  model  with  assump- 
tions, definition  of  parameters  and  system  variables  along 
with  some  basic  definitions  in  the  queueing  theory  forms  the 
first  half  of  this  chapter.   A  brief  discussion  of  the  pre- 
vious studies  on  this  model  then  follows. 

A.   DESCRIPTION  OF  THE  SYSTEM 

The  M/E, /c  queue  is  a  multi— server  queueing  system  with 
c  servers,  where  arrivals  occur  according  to  a  Poisson  pro- 
cess with  rate  A,  and  service  times  have  Erlang  distribution 
with  shape  parameter  k.   It  is  also  an  infinite  queue,  i.e. 
there  is  no  limitation  for  the  number  of  customers  in  the 
system. 

The  distribution  of  interarrival  times  X  is  given  by 

f  (x)  =  Xe~Ax  ,  x  >  0 

and  the  service  times  have  the  following  probability  density 
function: 

V*>  =  — H   yk_1e"y/6  ,    y  >  o 

(k-l)!3 


2 

with  mean  k3  and  variance  k$  . 
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1.   Definitions  of  Parameters 

a.  Arrival  Rate  and  Service  Rate 

Arrival  rate  to  the  system  is  the  reciprocal  of 
mean  interarrival  time  and  denoted  by  A.   Similarly,  service 
rate  of  a  server  is  defined  as  reciprocal  of  mean  service 
time  and  represented  by  y.   Then,  y  =  r-g   . 

b.  Traffic  Intensity 

Traffic  intensity  p  is  defined  by  the  following 
expression: 

p  =  A/cy  . 

c.  Offered  Load 

Offered  load  is  the  ratio  of  arrival  rate  A  to 
service  rate  y  and  denoted  by  a.  It  can  also  be  expressed 
by  product  of  number  of  servers  c  and  traffic  intensity  p 
as  follows 

a  =  cp  . 

d.  Coefficient  of  Variation 

The  ratio  of  standard  deviation  to  mean  is  de- 
fined as  coefficient  of  variation  of  a  probability  distribu- 
tion.  It  is  denoted  by  V.   For  Erlang  distribution  coeffi- 
cient of  variation  is  expressed  as 

v  ■  {k£?         1 
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2.   System  Variables 

a.  Number  in  the  System 

The  total  of  the  .number  of  customers  waiting  in 
the  queue  and  those  in  the  service  is  defined  as  number  in 
the  system  and  denoted  by  N.   The  number  of  customers  only 
in  the  queue  is  represented  by  Nq.   Expected  values  of  N 
and  Nq  are  denoted  by  L  (average  number  in  the  system)  and 
Lq  (average  number  in  the  queue) ,  respectively. 

b.  Waiting  Time 

The  time  spent  in  the  queue  by  a  customer  is  de- 
fined to  be  the  waiting  time  of  a  customer.   The  time  spent 
in  the  system  is  then  the  sum  of  waiting  time  and  service 
time  of  a  customer.   Tq  denotes  waiting  time  and  T  denotes 
the  time  spent  in  the  system.   Expected  values  of  T  and  Tq 
are  denoted  by  W  (average  time  in  the  system)  and  Wq  (aver- 
age waiting  time)  respectively. 

c.  State  Probabilities 

The  number  in  the  system  at  a  particular  time 
is  defined  as  the  state  of  the  system.   Then  state  probabili- 
ty P(N  =  n|t)  is  the  probability  of  the  system  being  in 
state  n  at  time  t. 

d.  Delay  Probability 

Delay  probability  is  the  probability  that  the 
arriving  customer  finds  all  the  servers  busy  and  enters  the 
waiting  line.   It  is  denoted  by  P (Tq  >  0).   Also,  by  defini- 
tion it  follows  that 

P(Tq  >  0)  =  P(N  ■     c)  . 
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3.  The  Steady— State  Solution 

The  steady— state  solution  is  defined  to  be  the  proba- 
bility distribution  of  the  number  in  the  system  when  the 
system  achieves  statistical  equilibrium.   In  the  steady— state, 
the  state  probabilities  do  not  depend  on  t,  i.e.  P(N  =  n|t)  = 
P (N  =  n) .   Also  it  follows  that 

00 

Z-   P(N  =  n)  =  1  . 
n=0 

In  queueing  theory,  it  is  a  well— known  fact  that  [1] 
the  steady— state  is  achieved  if  and  only  if  the  traffic  inten- 
sity p  is  less  than  1. 

4.  Assumptions 

The  M/E,/c  queueing  system  under  study  is  assumed  to 
have  the  following  properties: 

(i)    There  is  only  one  waiting  line  no  matter  how 
many  servers  there  are. 

(ii)   The  customer  at  the  head  of  the  waiting  line 
starts  getting  service  immediately  whenever  a  server  becomes 
available. 

(iii)  Order  of  service  is  first— come  first— served. 

(iv)   If  an  arriving  customer  finds  all  servers  busy, 
he  joins  the  waiting  line.   The  waiting  line  has  unlimited 
capacity. 

(v)    The  service  channels  (servers)  are  indexed 
consecutively.   If  an  arriving  customer  finds  more  than  one 
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server  vacant,  he  goes  to  the  Server  with  smallest  index 
(this  assumption  doesn't  cause  any  loss  of  generality,  but 
is  convenient  for  the  computer  simulation  model) . 

(vi)   The  servers  are  homogeneous,  i.e.  the  service 
distribution  is  the  same  for  all  servers. 

(vii)  Only  one  customer  at  a  time  can  be  served  by 
a  server. 

B.   FORMER  STUDIES 

There  are  numerous  studies  in  the  literature  of  queueing 
theory,  done  for  the  steady—state  solution  of  M/E,/c  queue. 
But  most  of  them  cover  only  some  particular  values  of  k  and 
c,  and  their  results  cannot  be  generalized.   Therefore,  such 
studies  will  not  be  mentioned  here  individually. 

The  earliest  suggestion  known  by  the  author  about  the 
solution  of  general  M/E, /c  system  was  mentioned  by  Lee  [2] . 
He  stated  referring  to  a  case  study  done  in  1956  that  mean 
waiting  time  can  be  approximated  by  use  of  David  Owen's  sug- 
gestion.  The  formula  given  in  [2]  was  said  to  be  usable  for 
0.7  _<  V  <_  1.0  and  as  follows, 

Wq  =  I  Wql  (1+V2) 

where  W   :  mean  waiting  time  for  M/E, /c  queue 

W  , :  mean  waiting  time  for  M/M/c  queue 

V   :  coefficient  of  variation  of  Erlang  (k) 

2 

distribution,  V  =  1/k. 

The  validity  of  the  approximation  was  demonstrated  by 
comparing  its  results  with  simulation  results  for  different 
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cases.   However,  the  approximation  is  just  for  mean  waiting 
time  or  queue  length,  hence  it  isn't  possible  to  approximate 
state  probabilities.   Also,  no  reference  was  given  for  de- 
tails of  David  Owen's  suggestion.   This  approximation  will 
be  mentioned  again  later  in  this  study. 

The  steady— state  solution  of  the  general  M/E,/c  queue 
was  first  studied  by  Mayhugh  and  McCormick  [3] .   The  results 
can  be  applied  to  the  cases  with  any  value  of  k  and  c.   How- 
ever, the  computation  procedure  is  so  complex  that  it 
requires  a  very  considerable  amount  of  work. 

A  parallel  study  was  also  done  by  Heffer  [4] .   The  solu- 
tion method  proposed  by  Heffer  differs  from  Mayhugh  and 
Mccormick's,  but  it  also  is  very  complex. 

The  reader  is  referred  to  the  references  [5]  and  [6]  for 
more  detailed  discussion  of  these  two  studies. 

The  most  recent  study  was  done  by  Yu  [5] .   It  concerned 
the  E  /E,/c  queue  with  heterogeneous  servers,  but  results 
were  also  stated  for  homogeneous  server  case.   Then  letting 
m  =  1  gives  the  solution  procedure  for  M/E, /c  queue.   However, 
like  the  other  two  studies,  the  computations  required  still 
demand  an  enormous  amount  of  work. 

The  theoretical  results  obtained  by  Heffer  [4]  and  Yu  [5] 
formed  the  basis  of  the  only  numerical  tables  available  for 
M/E, /c  queue,  prepared  by  Hillier  and  Lo  [6] .   The  tables 
have  cases  of  M/E, /c  queue  with  limited  values  of  k  and  c 
(k  <  8,  c  <  10),  and  a  few  cases  for  E  /E,/c  queue.   These 
tables  will  be  discussed  in  detail  in  the  following  chapter. 
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III.   DISCUSSION  OF  NUMERICAL  TABLES  AND 
DEVELOPING  A  NEW  COMPUTATION  METHOD 

In  this  chapter,  a  detailed  review  and  discussion  of  the 
numerical  tables  given  in  [6]  will  be  made.   Then,  using  the 
results  of  these  discussions,  a  new  computation  method  v/ill 
be  developed. 

Sample  pages    from  the   numerical  tables  are  given  in  Appendix  B. 

A.   DISCUSSION  OF  NUMERICAL  TABLES 
1.   General  Description 

The  tables  under  consideration  were  designed  for  gen- 
eral E  /E, /c  systems.   The  analysis  however  will  be  focused 
on  the  cases  in  which  m  =  1  (i.e.  M/E,/c  system) ,  which  forms 
a  large  porportion  of  the  tables.   As  mentioned  earlier,,  the 
values  of  k  and  c  for  various  tables  are  as  follows: 

k=2;  c=l,2,...,  10 

•K — -j  J      C  — "X  ;  —  •••♦/    D 

k=4;  c=2,3 

k=5,6, ... ,8;  c=2 
The  information  given  in  each  table  is: 
(i)    State  probabilities,  P(N=n),  and 
(ii)   Cumulative  state  probabilities,  P (N<n) ,  for 

il   -i-  f   <-  f   •  •  *  *r 

(iii)  Delay  probability,  P(T  >0); 

Si 

(iv)   Expected  number  in  the  system,  L; 
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(v)    Average  queue  length,  Lq; 

(vi)   Average  queue  length  for  the  equivalent  M/M/c 
system  with  the  same  arrival  and  service  rates,  L  ,  . 

(vii)  The  ratio  L  /L  , . 

All  this  information  is  given  in  each  case  and  for 
the  following  values  of  traffic  intensity  p  : 

p  =  0.10,  0.20,...,  0.50,  0.55,...,  0.95,  0.98,  0.99. 

So,  all  the  information  is  given  for  steady— state  con- 
ditions. 

Average  queue  length  L  ,  for  the  equivalent  M/M/c  sys- 
tem is  computed  by  using  the  following  formula  [7] , 

L  ,  =    (Cp)C?  Pn  (3.1) 

91   c!(l-p)2  ° 


where  c— 1      . 

V   (cp)3     (cp)c 


P 


o    j  =  0 


j!    +  cl(l-p) 


The  rest  of  the  information  given  in  the  tables  was 
computed  by  the  results  of  theoretical  work  mentioned  in  the 
last  chapter,  done  by  Heffer  and  Yu. 
2.   Use  of  Tables 

The  kind  of  information  listed  above  can  be  obtained 
directly  from  the  tables.   However,  by  making  use  of  some 
general  relationships  in  the  queueing  theory,  some  other  in- 
formation can  be  obtained  besides  those  mentioned  above. 

a.   Interpolation 

No  interpolation  method  is  specified  in  the  ex- 
planation sections  for  the  tables  in  [6] ,  for  the  values  of 
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p  different  from  those  given  above.   The  default  position 

taken  is  linear  interpolation. 

b.   Waiting  Times 

The  average  waiting  time  W  and  average  time  in 

Si 

the  system  W  are  frequently  of  interest  in  the  analysis  of 
queueing  systems.   These  values  can  be  obtained  from  the 
tabulated  values  of  L  and  L  respectively,  by  using  the  fol- 
lowing relationships  [8] , 

L 

W  =  -2.  w  =  - 

Wq    X     '  W    X  * 


The  probability  P(T  >t)  that  waiting  time  exceed- 
ing t  is  also  important  in  the  analysis  of  queueing  systems. 
Although  this  kind  of  information  is  not  given  in  the  tables, 
P(T  >t)  can  be  approximated  from  the  state  probabilities, 

Si 

P(N=n),  as  follows: 


P(T  >t)  =  Z^    P(N=n)  P{D< [k(n-c+l)-l] },  t>0, 
"  n=c 


(3.2) 


where  D  is  a  Poisson  random  variable  with  mean  ckyt.   The 
reader  is  referred  to  the  discussion  in  [6]  for  the  stochas- 
tic reasoning  of  this  relationship.   The  quantity 
P{D<[k (n— c+1)  — 1] }  can  be  obtained  from  Poisson  distribution 
tables  with  mean  ckyt.   If  this  mean  exceeds  25,  then  the 
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normal  distribution  with  the  same  mean  and  variance  can  be 
used  as  an  approximation  to  Poisson  distribution. 
3.   Discussion 

a.   Delay  Probabilities 

Similar  to  M/E,  /c,  M/M/c  denotes  a  multi— server 
queueing  system  with  Poisson  arrivals  of  rate  X    and  exponen- 
tial service  times  of  service  rate  y.   The  delay  probabili- 
ties for  M/M/c  queue  can  be  computed  by  using  the  formula 
[7], 


p(Tq>0)  =  ^TTI^T  Je=I = (3-3) 


(  y  (cP)j     (cP)c  \ 


A  chart  is  also  available  in  Appendix  A  for  c=2, 
3,... ,15  and  0<p<l,  for  P(T  >0). 

A  comparison  of  delay  probabilities  of  M/E, /c 
and  M/M/c  systems  for  selected  values  of  p,  c  and  k  is  given 
in  Table  I  of  this  study.   A  careful  examination  of  this 
table  shows  that  the  corresponding  delay  probabilities  of 
the  two  systems  for  the  same  value  of  p  are  very  close. 
Since  the  formula  or  charts  are  available  for  M/M/c  system, 
this  comparison  shows  that  they  can  also  be  used  to  estimate 
delay  probabilities  for  M/E, /c  system  with  a  very  small 
error,  usually  less  than  0.01. 
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b.   Average  Queue  Length  and  Number  in  the  System 
A  careful  examination  of  the  "Ratio"  column  of 
Table  I  in  [6]  indicates  that 


- 

(1  +  ±)  =  ±(1  +  V')        (3.4) 


lim  ^q  _  k+1    ln    1.    1,, 
p+1  L   "  2k~  "  ' 


which  is  essentially  the  same  as  Owen's  suggestion  [2],   How- 
ever, it  works  beyond  the  limits  of  coefficient  of  variation 
V  that  were  stated  by  Owen.   Also  (3.4)  is  true  even  when  p 
is  markedly  less  than  one.   The  theoretical  verification  for 
(3.4)  was  developed  by  Yu  in  [9]. 

A  more  precise  approximation  formula  was  also 
developed  by  Hillier  and  Lo  in  [6]  which  provides  more  accu- 
racy for  small  values  of  p.   In  this  case,  the  approximation 
formula  is  given  by 


Lq  *  |(l+g)  (1  +  |)  Lql  (3.5) 


where  g  =  f(p,k,c).   Then  the  expression  for  g  has  been  ob- 
tained by  linear  regression.   The  exact  coefficients  for 

g(p,k,c)  are  given  in  [6],   However,  a  rougher  expression 

«. 
which  is  more  convenient  for  computations  will  be  used  here, 

given   as 


g  =  xj(£+t>     (c-D2/3    {(1-p)    +    (1-p)2}       (3.6) 
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The  main  purpose  of  Hillier  and  Lo  for  introduc- 
ing g  into  approximation  formula  (3.5)  was  to  provide  means 
of  approximation  for  the  cases  with  larger  k  and  c  values 
which  were  not  covered  in  [6] ,  namely,  extrapolation  for  lar- 
ger values  of  c  and  k.   One  way  to  check  the  validity  of  the 
computational  formula  for  g  is  to  first  let  k  go  to  infinity, 
then  compare  the  values  of  average  waiting  time  W  computed 
by  the  formula  given  below  with  average  waiting  time  W   of 
M/D/c  queue.   This  latter  queue  has  Poisson  arrivals  with 
rate  X  and  constant  service  times  of  length  1/y  (or  k$) ,  and 
one  can  use  the  fact  that  the  distribution  of  service  times 
can  be  approximated  by  constant  service  time  distribution 
for  very  large  k  as  mentioned  in  Chapter  I. 

The  computation  formula  for  W  _  of  M/D/c  queue 
is  given  in  [2]  as  follows, 


=  I 

i=l 


W        =     L< 
qD 


— ia 


oo  oo 


(ia) 


:=ic 


j=ic+l      3 


i  l 


However,  a  more  convenient  form  for  computations  can  be  ob- 
tained as 


w 


qD 


i=l 


— ia , .  x ic 
e    (ia) 


TTcTT 


+  (1 


-i,  I 


— ia  ,  .  x 
e    (ia) 


p'j=ic+l     ^ I 


(3.7) 
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which  permits  the  use  of  Poisson  distribution  tables.  But 
it  still  requies  a  great  deal  of  computational  work.  For- 
tunately, a  very  convenient  chart  was  developed  by  Shelton 

[10]  for  W  _  for  l<c<100  .and  0.10<p<0.96. 

J       qD      —  — 

The  comparison  mentioned  above  can  be  made  hav- 
ing the  necessary  tools  available  and  Table  II  can  be  used 
for  this  purpose.   The  values  of  W   for  M/D/c  system  were 
obtained  from  Shelton 's  charts.   Average  waiting  time  W   for 
M/E,  c  system  is  computed  by  using  the  following  formula, 


wq  ■  ^t1  +  Ef(Eir)t^1)2/3{(1-<»  +  <i-p)2}]<i+  iKi 

where  L  ,  is  computed  from  (3.1).   If  the  expression  given 
by  (3.6)  is  true  for  very  large  values  of  k,  then  it  should 
also  be  true  that 


lim  Wq  =  "qD 


Taking  the  limit  of  W  gives 

lim  Wg  =  ^[l  +  ^  (c-l)2/3  {(1-p)  +  (1-p)2}]  Lgl 
k+oo 

On  the  other  hand, if  g  is  omitted,  the  limit 

will  then  be 

lim  W  =  tt  L  ,  =  -7T-   W  1 
q    2A   ql    2   ql 

k-K» 
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A  comparison  of  W   and  these  two  limits  for 
eight  different  cases  is  given  in  Table  II.   Investigation 


TABLE  II. 

> 

Comparison 

of  W  .  and 
qD 

Limits  of 

W  With  or 

q 

Withou 

Case 

:  No 

c 

P 

2   ql 

W  (k=°°) 

q 

V 

1 

2 

0.80 

1.185 

1.209 

1.08 

2 

2 

0.80 

1.422 

1.451 

1.296 

3 

5 

0.90 

2.287 

2.34 

2.16 

4 

5 

0.90 

1.144 

1.17 

1.08 

5 

10 

0.90 

1.337 

1.391 

1.36 

6 

10 

0.90 

1.505 

1.563 

1.53 

7 

20 

0.90 

3.305 

3.766 

3.36 

8 

20 

0.90 

5.508 

5.88 

5.6 

of  these  results  shows  that,  for  almost  all  the  cases  the 

limit  with  g  omitted  is  closer  to  W  _  than  the  other  limit. 

qD 

This  indicates  that  the  formula  (3.6)  should  be  reviewed 
for  large  k  values.   However,  this  revision  can  be  done  only 
after  some  additional  tables  become  available  for  the  cases 

not  covered  in  [6] . 

«. 
c.   Waiting  Times 

The  method  of  approximating  the  probability 

P(T  >t)  that  the  waiting  time  will  exceed  t  using  the  state 

probabilities  P (N=n)  was  described  earlier  in  the  chapter. 

However,  this  method  assumes  that  an  arriving  customer  will 
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have  to  wait  for  approximately  k(n— c+1)  service  phase  (ex- 
ponentially distributed  with  mean  $)  completions  before  a 
server  becomes  available  to  start  his  service.   It  is  stated 
in  [6]  that  this  approximation  for  P(T  >t)  should  be  better 
for  larger  values  of  p  and  t  due  to  this  assumption.   It 
should  be  added  that,  for  large  k  and  small  t,  the  approxi- 
mation fails  greatly  in  representing  the  actual  distribution 
of  waiting  times  for  the  same  reason. 

B.   NEW  COMPUTATION  METHOD 

The  construction  of  a  practical  computation  method  using 
the  facts  obtained  from  the  preceding  discussions  will  form 
the  rest  of  the  chapter.   This  new  method  is  desired  to  be 
applicable  to  the  queueing  problems  involving  M/E,/c  systems 
met  in  practice,  with  simple  and  quick  computational  work 
without  any  need  for  tables,  interpolations  etc. 

The  kinds  of  information  which  are  desired  to  be  able  to 
compute  by  the  new  method  for  M/E, /c  queueing  system  are 
mainly 

(i)    Delay  probability, 

(ii)   Average  waiting  time,  average  time  in  the  system, 
average  queue  length  and  average  number  in  the  system, 

(iii)  Probability  distribution  of  delay  times  and  proba- 
bility that  waiting  time  will  exceed  a  specified  time  length, 

1.   Delay  Probability 

In  the  discussion  earlier  it  was  pointed  out  that 

the  delay  probabilities  P(T  >0)  for  the  systems  M/Ev/c  and 

q  K 
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M/M/c  are  very  close  for  the  same  value  of  p.   An  examina- 
tion of  the  comparisons  given  in  Table  I  shows  that  the 
delay  probabilities  are  very  close  for  M/M/c  (k=l) ,  M/E, /c 
(k>l)  and  M/D/c  (k=°°)  systems  with  the  same  value  of  p. 
The  differences  between  the  delay  probabilities  for  the  same 
p  are  less  than  0.01  for  almost  all  the  cases  and  gets  even 
smaller  for  p  close  to  one.   This  comparison  suggests  that 
delay  probabilities  of  the  M/M/c  system  can  be  used  for  cor- 
responding M/E,/c  system  also.   These  delay  probabilities 
are  computed  by  using  formula  (3.3).   The  formula  is  not 
suitable  for  large  values  of  c,  however,  and  a  chart  is  given 
in  Appendix  A  for  up  to  about  16.   For  greater  values  of  c, 
the  charts  are  also  available  in  [10]  and  [11] . 

2.   Average  Queue  Length  and  Average  Waiting  Time 

In  the  previous  discussion,  after  examination  of  the 
tables  in  [6]  it  was  concluded  that  the  average  queue  length 
for  M/E, /c  system  can  be  approximated  by 

Lq  =  i(l+g)  (1  +  l)Lql  (3.8) 


where  L  ,  is  computed  according  to  (3.1),  and  g  is  computed 
by 

g  =  TZ    kTT(c"1)2/3  {(1""P}  +  <!-P)2>  • 

Contribution  of  g  in  equation  (3.8)  is  negligible 
for  practical  purposes  in  most  cases.   Also  it  was  shown 
earlier  in  this  chapter  that  the  formula  g  requires  some 
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modifications.   Therefore  it  will  be  omitted  in  further 
discussions. 

Now,  attention  will  be  focused  on  equations  (3.1) 
and  (3.3).  Rewriting  equation  (3.1)  in  a  different  form, 
then  substituting  (3.3)  gives 


(cp) 


ql    (1-P) 


c!  (1-p) 


rc-1  i       c  " 

£  <cp>  +  (cp) 

j=0  1[  c!  (1-p) 


and 


p«P(T  >0) 

~d 


ql     (1-P 


Then  it  follows  that 


L  ,    P  T  >0 

w      =  _ai  =   <3 

ql    X    cy(l-p) 


Using  (3.8)  with  g=0,  some  approximation  formulas 
for  average  queue  length  and  average  waiting  time  are  ob- 
tained respectively  as  follows 


,         P(T   >0)  , 

Lq  -  h  p  -tSpt  (1  +  e 
i     pcyo)         i 

Wq    =   7         cTrt=pT    U    +    k> 


(3.9) 
(3.10) 


where  P(T  >0)  is  delay  probability  for  M/M/c  system  which 
can  be  obtained  from  the  charts  very  easily.   Then  it  is 
very  simple  to  compute  L   and  W  using  equations  (3.9)  and 
(3.10)  . 
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The  average  number  in  the  system  L  and  average  time 
in  the  system  W  can  be  computed  by  following  well— known  rela- 
tionships of  queueing  theory, 

*  -  Lq  +  "J"  '  "  "  X  ' 

3.   Waiting  Time  Distribution 

Suppose  that  the  conditional  waiting  time  distribu- 
tion given  T  >0  for  M/E,/c  queueing  system  can  be  approxi— 
q  k 

mated  by 


•L-  .1. 

P(T  >t|T  >0)  =  e 


Then  using  Bayes'  theorem  it  follows  that 


P(T  >t)  =  P(T  >t|T  >0)P(T  >0)  =  e  bt  P  (T  >0)  .  (3.11) 
h        h    ^d      y  q 


Let  Fm  be  the  CDF  of  waiting  times. 
Tq  3 


Then 


PTq(t)  =  1-P(T  >0)  e  bt  ,         tO 


Now  average  waiting  time  can  be  computed  as  follows 

00  00 

Wq  =  E(Tq)  =  Qf    (l-FTg(t))dt  =Q/  P(Tq>0)e"btdt 

P(T  >0) 
Wq  =  1 .  (3.12) 

Suppose  the  parameter  b  is  modeled  as 

b  m   2cy (1-p) 
(1+V2) 
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Then, 

P(T   >0)  0 

Wq   =   2cu?I=pJ U+V    '  (3-13) 


where  V  is  coefficient  of  variation  of  service  time  distri- 
bution. Notice  that  equation  (3.13)  is  exactly  the  same  as 
equation  (3.10).  However,  there  can  exist  some  other  wait- 
ing time  distribution  to  give  the  same  average  waiting  time 
as  given  by  (3.13),  so  it  must  be  shown  that  for  M/E,/c  sys- 
tem, the  distribution  of  waiting  times  can  be  approximated 
by 

-  2cy (l-p)t 
F   (t)  =  1-P(T  >0)e    (1+V2)     ,  t>0       (3.14) 

If  (3.14)  is  true,  then  P(T  >t)  computed  by  (3.11) 
should  be  approximately  equal  to  the  one  obtained  by  the  pro- 
cedure suggested  by  Hillier  and  Lo  [6]  which  uses  state  pro- 
babilities as  described  earlier  by  equation  (3.2).   Table 
III  gives  the  comparison  of  P(T  >t)  values  obtained  by  these 
two  different  formulas  for  various  values  of  t.   It  should 
be  recalled  that  approximation  by  state  probabilities  is 
poor  for  small  values  of  t.   However,  for  t>W  the  two  re- 
sults agree  very  well.   The  difference  is  ease  of  computa- 
tions.  Equation  (3.11)  is  much  easier  than  equation  (3.2)  - 
to  compute  the  same  value. 

The  validity  of  approximation  (3.14)  for  small  values 
of  t  will  be  shown  by  simulation,  which  forms  the  next  chap- 
ter.  Then  the  suggested  computation  method  will  formally 
be  described  on  a  step— by— step  basis  in  last  chapter. 
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IV.   COMPUTER  SIMULATION 

A  computer  simulation  model  for  M/E, /c  system  is  used 
to  investigate  the  validity  of  approximation  for  waiting 
time  distribution  introduced  in  the  preceding  chapter. 
Validity  of  this  approximation  for  moderate  and  large 
values  of  t  was  demonstrated  in  the  same  chapter.   There- 
fore, analysis  will  now  be  focused  on  the  approximation 
for  small  values  of  t. 

A.   COMPUTER  PROGRAM 
1.   Model 

Since  the  model  used  to  construct  the  simulation 
program  is  based  on  the  same  assumptions  as  stated  in  Chap- 
ter II,  it  won't  be  described  again  in  this  chapter.   The 
computer  language  selected  was  SIMSCRIPT  which  is  especial- 
ly convenient  for  simulation  of  multi— server  queueing  sys- 
tems.  The  reader  is  referred  to  [12]  and  [13]  for  a  brief 
idea  about  SIMSCRIPT  if  he  is  not  familiar  with  this 
language.   Main  reference  for  SIMSCRIPT  however  is  [14] . 

One  way  to  check  for  small  t  the  validity  of  approxi- 
mation distribution  developed  earlier  is  to  keep  the  fre- 
quency table  of  waiting  times  of  the  customers  during  the 
simulation.   Let  M(t)  be  the  number  of  customers  with  wait- 
ing time  greater  than  t  and  M  be  total  number  of  customers 
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served  during  simulation  period.   Then  the  estimator  for 
probability  of  waiting  time  greater  than  t  is 


P(Tq>t)  =  ^1  (4.1) 


This  value  can  be  compared  with  the  value  given  by  the  ap- 
proximation formula  to  check  its  validity. 

The  computer  program  is  given  in  Appendix  B.   It  is 
written  in  the  way  that  it  would  be  possible  to  obtain  the 
estimates  of  all  kinds  of  information  about  the  queueing 
system  being  simulated;  average  waiting  time  and  average 
queue  length,  state  probabilities,  delay  probability  etc. 
However,  only  the  estimates  given  by  (4.1)  will  be  discussed 
here. 

2.   Variance  Reduction 

One  of  the  problems  encountered  in  the  simulation 
of  queueing  systems  is  the  high  variability  of  waiting  times, 
especially  when  traffic  intensity  p  is  high.   Also  strong 
positive  correlation  between  the  waiting  times  of  consecu- 
tive customers  causes  serious  errors  if  the  standard  for- 
mula is  used  to  estimate  the  sample  variance  of  waiting 
times  since  this  formula  will  underestimate  the  true  vari- 
ance [1]  . 

There  are  several  methods  developed  to  overcome  this 
difficulty  with  high  variability.  Interested  reader  can  find 
comprehensive  information  in  references  [12] ,  [15]  and  [16] . 
The  variance  reduction  method  which  will  be  used  here  is 
antithetic  variates. 
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In  the  simulation  program,  two  different  random  num- 
ber streams  are  used  to  generate  interarrival  and  service 
times.   Once  a  uniform  variate  between  0  and  1  is  generated, 
then  the  random  variate  from  a  desired  distribution  can  be 
obtained  by 

X  =  F "X(u)  (4.2) 

where  u  is  a  uniform  variate  between  0  and  1  and  F    denotes 

x 

the  inverse  of  CDF  of  X  [12] .   Then,  for  example,  a  sequence 
of  exponentially  distributed  random  variates  with  mean  1/X 
can  be  obtained  from  a  stream  of  uniform  variates  by  using 

X  m   ZJL  in(u)  (4.3) 

Let  Z,  and  Z2  be  estimators  of  a  parameter,  obtained 
from  two  different  simulation  samples  and  Z3  =  -j  (Z,  +  Z2) . 
Then 

Var(Z3)  =ivar(Z1+Z2)  =i  [Var  (Z^  +Var  (Z2)  +Cov  (Z^  Z2)  ] 

(4.4) 
which  implies  that  maximum  negative  correlation  between  Z.. 
and  Z2  would  minimize  Var(Z3).   This  can  be  obtained  by  us- 
ing 1— u  in  (4.2)  in  place  of  u  for  random  variate  genera- 
tion in  the  second  simulation  run. 

During  the  simulation,  a  random  sequence  of  large 
service  times  or  short  interarrival  times  would  cause  long 
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waiting  times  and  conversely.   Using  the  procedure  described 
above,  a  sequence  of  large  waiting  times  would  become  sequence 
of  short  waiting  times  in  the  second  run;  in  other  words,  the 
waiting  times  in  two  different  runs  will  be  negatively  corre- 
lated.  Let  P. (t)  be  the  probability  that  waiting  times  will 
exceed  t  in  the  ith  run,  i  =  1,2.   Using  antithetic  variates, 
P, (t)  and  P2 (t)  will  be  negatively  correlated?  then  the  esti- 
mator will  be  computed  by 


P(T  >t)  =  |(P1(t)  +  P2(t)) 


where  P,  (t)  and  P?^)  are  computed  according  to  (4.1). 

Generation  of  antithetic  exponential  variates  for  in— 
terarrival  times  is  given  by  (4.3).   However,  it  is  not 
straightforward  to  generate  antithetic  Erlang  variates  since 
the  CDF  cannot  be  inverted  in  closed  form.   Generating  Erlang 
variates  as  sum  of  k  exponential  variates  with  mean  (3  does 
not  help  since  equation  (4.2)  cannot  be  utilized  to  obtain 
antithetic  variates.  Nevertheless,  one  way  to  solve  this 
problem  is  to  store  the  CDF  table  of  chi— square  distribution 
in  an  array,  then  compute  F  ""  (u)  by  linear  interpolation  to 
get  a  chi— square  random  variate,  and  obtain  the  Erlang 
variate  by  the  following  transformation 

Ek  =  2  BCk 
where  E,  is  the  Erland  variate  with  mean  k$  and  C,  is  the 
chi— square  variate  with  degrees  of  freedom  k.   The  tables 
for  CDF  of  chi— square  distribution  are  available  in  [17] . 
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3.   Data  Collection 

a.   Selection  of  Input  Parameters 

The  input  values  for  each  pair  of  runs  were  arbi- 
trarily selected  to  give  certain  traffic  intensity  p  values, 
usually  0.8  and  0.9  in  order  to  get  a  wider  range  of  waiting 
times  even  though  they  cause  high  variability.   The  width  of 
frequency  intervals  was  chosen  to  be  0.25  min.  so  that  f , , 
for  example,  would  show  the  number  of  customers  with  waiting 
time  less  than  0.25  minutes.   Then  M(t) ,  the  number  of 
customers  with  waiting  time  greater  than  t  can  be  computed 
by 


:t)  =  S 


M(t)  =  Z-r     f, 
i=4t+l   1 


since  t  will  be  in  multiples  of  0.25.   This  computation  is 
done  in  the  computer  program. 

b.   Initial  and  Final  Conditions  of  the  System 

Initial  and  final  conditions  of  the  simulated 
system  are  important  since  they  effect  the  value  of  para- 
meter estimated.   The  approach  taken  in  this  study  was  to 
use  epochs.   If  the  time  at  which  the  number  in  the  system 
N  changes  from  0  to  1  by  an  arrival  is  defined  as  a  regenera- 
tion point,  then  the  interval  between  two  successive  regene- 
ration points  can  be  defined  as  an  epoch  [12] .   It  is 
obvious  that  the  sequences  of  waiting  times  in  two  different 
epochs  will  be  independent  from  each  other.   The  epochs  tend 
to  be  lengthy  with  high  traffic  intensity  p,  large  arrival 
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rate  X,  with  large  number  of  servers  c,  or  a  combination  of 
these  three  factors.   Before  each  simulation  run,  the  number 
of  epochs  for  which  the  simulation  is  to  be  run  was  deter- 
mined as  an  input  considering  those  three  factors.   Exper- 
ience showed  that  the  first  few  epochs  tend  to  have  waiting 
times  smaller  than  usual,  so  they  were  omitted  for  data 
analysis.   This  isn't  feasible  for  cases  with  long  epochs; 
however  the  first  couple  of  epochs  have  enough  information. 
The  simulation  runs  were  started  with  N  =  0  and  an  arrival 
so  that  starting  time  would  be  a  regeneration  point  and 
stopped  after  the  number  of  epochs  determined  earlier  is 
completed,  i.e.  N  =  0  was  the  final  condition. 

P, (t)  or  P2(t)  are  given  in  the  output  of  program 
for  t  =  0.25,  0.50,  0.75,  ....  etc.   Then  to  get  the  estima- 
tors P(T  >t) ,  all  one  has  to  do  is  to  average  them  for  cor- 
responding t  values. 

B.   RESULTS 

A  comparison  of  waiting  time  probabilities  obtained  from 
simulation  and  approximation  method  is  given  in  Table  IV. 
Investigation  of  this  table  indicates  that  the  approximation 
method  agrees  quite  well  with  the  results  of  simulation  for 
small  t  values.   The  difference  between  the  corresponding 
values  for  moderate  or  large  values  of  t  for  which  the 
validity  of  approximation  was  demonstrated  earlier  explains 
the  difference  between  the  values  when  t  is  small. 
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V.   CONCLUSION 

It  was  shown  in  the  last  chapter  that  the  approximation 
method  developed  earlier  for  the  distribution  of  waiting 
times  can  also  be  used  for  small  values  of  t.   The  method 
for  computing  steady— state  information  of  M/E,/c  queueing 
system  can  now  be  formally  described. 

A.   DESCRIPTION  OF  THEPROPOSED  COMPUTATION  METHOD 

After  analysis  of  data  and  the  decision  has  been  made 
that  the  particular  system  under  study  can  be  approximated 
by  M/Ek/c  system  as  described  in  Chapter  II,  and  also  having 
the  estimators  for  X,    k  and  3  obtained  (by  maximum  likeli- 
hood or  method  of  moments) ,  one  is  ready  for  computations 
to  get  the  desired  information  for  the  system. 
1.   Delay  Probability 

First  compute  service  rate  y  as 

1 

y  =  k3  * 

Then  compute  traffic  intensity  p  and  offered  load  a  by 

X  ,  X 

p  =  —         and       a  =  — 

cu  y 

respectively.   Then  enter  the  chart  in  Appendix  A  with  a  to 
get  delay  probability  P (T  >0) .   Use  the  charts  in  [10]  or 

Si 

[11]  if  number  of  servers  c>16. 
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2.   Average  Waiting  Time  and  Average  Queue  Length 
First  compute  a  system  constant,  namely  b,  by 

k  -  2cy(l-p) 
1 
1  +  K 


Then  average  waiting  time  W   is  computed  as  follows 

si 


P(T  >0) 

w  =   9 

q       b 


Also,  average  queue  length  L   is  given  by 

si 

L   =  AW   . 

q    q 

3.  Average  Time  in  System  and  Average  Number  in  System 
Average  time  in  system  W  and  average  number  in  sys- 
tem L  will  be  computed  by 

W  =  W  +  i       and      L  =  L  +  - 

q   y  q   y 

respectively. 

4.  Waiting  Time  Distribution 

The  CDF  of  waiting  times  is  given  by 

FTq(t)  =  l-P(Tq>0)e"bt  ,       t>0 

where  b  is  the  system  constant  computed  in  step  2. 

Then  the  following  probabilities  can  be  computed 

(i)    p(Tq>ti)     "  P(T  >0)e~btl 

(ii)   P(t,<T  <t  )  =  P(T  >0) (e~btl  -  e"bt2) 

"  "  — bt 

(iii)  P(T  <t,)     =  1-P(T  >0)e    1 

*i   •*•  SI 

as  needed  by  the  user. 
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B.   DISCUSSION 

1.  Advantages  of  the  Method 

(i)    Needless  to  say,  the  most  important  advantage 
of  the  method  is  simplicity.  The  kind  of  information  listed 
above  can  be  computed  in  minutes  given  c,  X,  k  and  \i   which 
would  be  needed  to  compute  anyway.   All  of  the  computations 
can  be  laid  down  on  one  regular  size  page  so  that  it  is 
very  easy  to  follow  by  somebody  else. 

(ii)   In  most  cases,  one  has  to  compute  the  above 
listed  information  to  determine  the  effect  of  changing  the 
number  of  servers  c  on  the  selected  system  variables.   This 
multiplies  the  savings  of  time  and  computational  effort. 

(iii)  No  tables  are  necessary  except  for  the  delay 
probability  chart.   No  interpolation  would  be  needed. 

2.  Disadvantages  of  the  Method 

(i)    Since  the  method  gives  an  approximation,  it 
is  not  too  precise  even  though  the  results  are  almost  always 
in  +  5  percent  of  the  actual  value,  and  in  +  2  percent  for 
most  cases. 

(ii)   The  approximation  of  state  probabilities  with 
this  method  does  not  appear  to  be  direct. 

3.  Applications 

One  of  the  characteristic  applications  of  the  queue— 
ing  theory  is  to  investigate  the  effect  of  changing  the 
parameters  or  number  of  servers  on  the  measure  of  effective- 
ness (MOE)  assigned  by  the  management.   This  measure  of 
effectiveness  can  be  delay  probability,  average  waiting  time 
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or  probability  of  waiting  time  exceeding  time  t  etc.   The 
marginal  utility  of  c  +  1st  server  will  be  the  difference 
between  the  values  of  MOE's  for  c  and  c+1  servers.   Then  the 
decision  criterion  whether  or  not  to  hire  the  c+lst  server 
is  how  his  cost  compares  to  his  marginal  utility. 

Another  managerial  objective  may  be  to  require,  for 
example,  the  probability  of  waiting  time  exceeds  4  min., 
P(T  >4)  =  0.05.   The  number  of  servers  necessary  to  achieve 
this  purpose  will  then  be  of  interest. 

It  is  possible  to  extend  the  variety  of  examples. 
The  common  property  of  the  above  examples  is  that  precision 
greater  than  that  of  suggested  approximation  is  not  needed 
to  make  the  decision. 
4.   Final  Remarks 

The  author  wishes  to  emphasize  that  without  the 
numerical  tables  provided  by  Hillier  and  Lo  [6] ,  it  would 
have  been  impossible  to  develop  such  an  approximation  method, 
The  method  may  need  some  modifications  as  new  tables  for 
the  cases  not  covered  in  [6]  become  available. 
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APPENDIX  A 

CHARTS  FOR  DELAY  PROBABILITY 

In  the  two  following  charts,  the  delay  probabilities 
can  be  found  for  M/M/c  system  which  were  shown  to  be  very 
close  to  those  of  M/E, /c  system,  with  number  of  servers  c 
up  to  about  16.   One  has  to  enter  to  chart  with  offered 
load  a  as  abscissa,  then  delay  probability  P (T  >0)  is  the 
ordinate  value  of  the  point  on  the  proper  c  curve  which 
has  a  as  abscissa  value. 
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APPENDIX   B 
SAMPLE    PAGES    FROM   NUMERICAL   TABLES 
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0. 8179*10 

0.  16*1) '9 

0.01464797 

0.0O1  U2!<I 

6.6841***-0* 
J.2**01T*-06 
l.600804*-07 
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3  I  *<"  I   I 

0. 1179210 
0.4820784 
0. 9467*64 
0.4444247 
0. 4444466 
0.9999998 
0.9*99999 


STATE 
I 

7  I 

8  I 

9  I 
10  I 
II 

12  I 

I) 


PIN-  I  I 

J.»44422*-04 
Z.ltT*»?*-l«i 
7.71*0  ■rf'-W 
2.7»4»28*-l3 
9.481  172'  -8* 
3.345)45' -16 
1.19027»*-17 

2  •  OxU 


»   1 

0.8652663 

1    1 

0.2698672 

1    1 

0.0*393797 

)   1 

0.0082)260 

1  1 

0.OOO98232O2 

)    1 

0.000102*607 

t    1 

4.794626'-06 

1   1 

8.S9080S,-07 

•  I 

7.83J906,-O8 

N      -       1             •   - 

inn 

i 

PIN-l  1 

0     1 

0.5332*89 

i    I 

0.3293022 

7     1 

0.1033992 

1      | 

0.02*288*6 

»     | 

0.00*626253 

5     I 

0.0007867696 

t     I 

O.OO01256372 

T      1 

l.93*6W,-0S 

1     1 

2.9*07*5'-0* 

I      •      I 

suit 

1 

eis-t  i 

0      1 

0.6233869 

1    1 

0.3332261 

i    1 

0.1353*5* 

)    1 

0.030058)2 

♦    1 

0.01)39216 

i    1 

0.003364750 

8      1 

0.0007926517 

7      1 

0.00018)1)28 

•      1 

6.290506'-05 

9      1 

9.6lB»29*-fl6 

18      1 

2.20Z*3S*-0* 

M       •       1             t 

mil 

i 

■MK»II 

0 

0.3265101 

1 

0.3*69791 

2 

0.19728*6 

3 

0.08362*06 

« 

0.030*2698 

3 

0.01023*98 

8 

0.0033)9760 

7 

0.00107)73* 

1 

|       0.009)6)92*2 

9 

|       0.00311009*0 

10 

1       ).S2*6+**-03 

11 

I       1.121*7*'-05 

H      •       1            • 

$t»r£ 

I 

HN-II 

0 

|              0.2829307 

1 

|              0.3)61386 

2 

|               0.21273*0 

3 

|               0.1018312 

8 

|            0.0*223055 

3 

I            0.016)*13« 

8 

I         0.0061*1*69 

T 

I         0.00228*550 

8 

I       0.0009*755*1 

9 

|       9.000)1*3163 

10 

I       0.0001165765 

11 

|       ».J2«122,-05 

R       •      1 

p  ( •»<  - 1  1 

0.6652663 
0.93*7336 
0.9906716 
0.99890*2 
0.9998868 
0.9999892 
0.9999990 
0.9999999 
0.9999999 
•      6 


STATE 
I 

9 
10 
11 
12 
1) 
16 
IS 
16 
17 
C 


PIN- I > 

|  *>.*0T7il« -04 

|  5.837*5)'-10 

|  3.038657--U 

J  *.6022l»,-l2 

|  3.112231*-13 

|  J.303994*-!* 

I  2.1*)789*-IS 

I  2.*428*0*-l6 
I       2.152598*-17 

•  2  •  BnO 


»<><<•;  i 

0.3332*89 
0.8667310 
0.9701303 
0.99**S8T 
0.9990650 
0.9998517 
0.9999772 
0.9999965 
0.9999995 


STATE 
I 


Pl>t»ll 

*.**8650*-97 
6.72*838'-38 
1.0170*0* -OS 
l.*389*)*-04 
2.329)48'-l0 
3. 526286* -11 
*.338370*-12 
8.0)16**'  -1  1 
1.223*53'-13 


PIN<-II 

0.9999999 
0.9999999 
0.4444444 
0.9499999 
0.9999999 
0.9999999 
0.9999999 


0.9999999 

0.9999999 

0.9949999 

0.  9999999 

0.999999) 

1.000030 

1.000000 

1.000130 

I. 000000 

.30 


0.94W444 
0.4944999 
9.9999999 
0.9999999 
0.9999999 
0.9999994 
0.9999999 
0.9999999 
0.9994999 


STATE 
I 

0  I 

1  I 

2  I 

3  I 


STATE 
I 


<MN<-I > 

0.*233869 
0.7766131 
0.9319585 
0.9820168 
0.9956089 
0.9989697 
0.999762* 
0.9999..55 
0.99*9875 
0.9999971 
0.9999993 
•       8  . 


C     •      I 


STATE 
I 


Rh9 


PIN'l ) 


PIN<«II 

0.3265101 
0.673*899 
0.87077*6 
0.95*3986 
0.93*8256 
0.4950306 
0.993*20* 
0.999*9*1 
0. 9998381 
0.9999*81 
0.999983* 
0.99949*6 


PIN<"I I 

0.2829307 
0.6170692 
0.8248033 
0.93165*6 
0.  9738651 
0.9902267 
0.99*3682 
0.9986528 
C.°995093 
0.  J9481** 
0.9994)12 
0.99997*4 


11  |  ».0*S1*3* 

12  I  I. 155876' 
1)  I  2.6*8*09* 

14  |  6.06932)* 

15  I  1.390**** 

16  I  3.183936* 
IT  I  7. 299952" 

18  I  1.6726*)' 

19  I  3.8)2525' 

20  I  8.781*72' 

21  I  2.012101' 
C  -  2 

STATE 

J  PIN- 


■4/7 
-07 
-3  8 
-09 
-09 
-13 
-11 

-11 
-U 
-13 

-IS 

MHO 


12  1  3.613901' 

13  I  X. 15726)' 
16  |  3.7C546*' 

15  I  1.196716' 

16  I  3.80018*' 
IT  I  1.216919 
18  I  3.894893 

14  I  i.2*r«as 

20  I  3.99606) 

21  I  1.2  796*5 

22  I  *.04776O 

23  |  1.312219 


•      2 


-06 
-06 
-37 
-0? 
-08 
-03 
-39 
i-09 
'-13 
'-10 
'-li 
'-11 
TWO 


STATE 
I 

12  I 

13  I 
1*  I 
13  I 

16  I 

17  I 
IS  I 
19  I 
20 

21 
22 
23 
C      • 


PIM-I   I 

1.603445*-0* 
3.4»44*0'-)6 
2.2071l»*-06 
S.l»722o'-37 
3.0)1022*  -37 
1.126S72'-i»7 
*.178»7»*-81 
1.5501  75* -U8 
S.7*0318*-J9 
2.133059*-09 
7.412502'-10 
2.9351W-19 
2  .  <U*0 


pi~<t! t 

0.9994998 
0.9949449 
0.499999  4 
0.9999999 
0.9949949 
0. 1449444 
0.4994499 
0.9999999 
0.4949999 
0.9999944 
0.4999494 
.50 


PIN<**II 

0.4444433 
3.494999* 
0.9944998 
0.9944944 
0.4449999 
0.9499999 
0.9999099 
0.9994999 
0.9944444 
0.4949999 
0.4944944 
0.9944999 


PINOII 

0.9949903 
0.9949463 
0.9949917 
0.999999  5 
0.4949491 
0.9949999 
0.9999444 
0.9449444 
0.9949999 
0.9999949 
0.9949994 
0.4994944 


STATE 
t 

0  I 

1  I 

2  I 

3  I 


PIN-I » 

0.23*2829 

O.  241*)*! 

0.227545) 

0.1 159851 

0.0715*6*9 

0.03559238 

0.01733*83 

0.008)877*1 

0.00*092620 

0.091937785 

0.0009*58*77 

0.000*569793 

0.0002207896 

X      -      1  • 


P(Nrl) 

0.1*89107 

0.2621716 

0.22*9721 

0.1*81*35 

0.987)9370 

0.0*97426* 

0.02673223 

0.01*5413* 

0.0079*9122 

0.006)30911 

0.002)39518 

0.001285522 

0.000  7093893 

0.0093813431 

M      •       I 


PIIOll 

0.133*713 

0.2286572 

0.21*25*7 

0.15621** 

0.1011311 

0.063**116 

0.03841735 

0.02377324 

0.01*510)3 

0.009956063 

0.005*05270 

0.00)29916) 

0. 09201 S69L 

0.001229035 

O.OOC750190) 

0.000*578896 


STATE 

I 


0.60 


STATE 
I 

0 
1 
2 
3 


PIN»I I 

9.2*22  59* 

9.)15*912 

0.22)17)6 

0. 1197532 

0.036125*1 

0.02*71790 

0.0106191) 

0.00*52*238 

O.OOl"  >409 

0.0098179**1 

0.000)-777*9 

0.0091*7977* 

6.288060'-0» 


PINOI  > 

0.2*2254* 
0.557  7*06 
0.79C4l»2 
0. 90C6675 
9.9587929 
0.9815107 
0.9921238 
0.49665)1 
0.  94HS7  70 
9.499)9*4 
0.9447*27 
0.4449906 
0.994953* 


STATE 
I 


PIN- I  I 


13  I  2.673e>) 

I*  |  I. 136478' 

15  I  *.8)*7J6 

16  I  2.0559)) 

17  I  8.7V.8JJ 
IS  I  3.717262 

19  I  1.580*69 

20  I  6.7<l!92 

21  |  2. 8*31. J 

22  I  1.213129 

23  I  5.1678/0 
2*  1  2.197502 
23  I  9.)**)05 


-05 
-05 
-0* 
-36 
-0  7 
-37 
-07 
-J8 
-J* 
-08 
>HJ9 
'  -0* 
'-10 


0.4994102 
9.9949915 
0.944996* 
9.9499995 
0. 0444943 
0.4444492 
0.9444494 
0.4944494 
3.4999999 
0. 9999999 
0.9994999 
9.9994444 
0.4999994 


P<NMI 

0.10*7049 

0.1903802 

0.14*3581 

0.1355139 

0.1120363 

0.077*9*80 

0.03299866 

0.03599627 

0.026*6579 

0.016632  78 

0.01139784 

0.997697837 

0.005226704 

0.003553*71 

0.002*15890 

0.0016*2*46 

0.031116673 

0.0CO7591899 

0.0005161*87 


PIN-I1 

0.075785*6 

0.1*8*290 

0.16*1668 

0.1*36229 

0.1138889 

0.08  706796 

0.06581*50 

0.0*461175 

0.9)7)8017 

0.92316)50 

0.0212147* 

0.91548012 

0.3120*636 

0.039076*13 

0.0068)8687 

0.905152658 

0.00  3392306 

0.032425151 

0.03220)476 

0.001960  601 

0.001241 192 

0.000**27142 

0.009719248) 

0.0015)5179* 

0.009-0)23*5 


STATE 

PIM<>  1 1 

1 

p  1  (••  1 1 

0.20*2829 

1)       1 

0.0001066  2*9 

0. '.457170 

1*       1 

S.l*»022*-3i 

0.723312* 

IS       1 

2.  *90l76'-05 

0.6542475 

16       1 

l.20)l33'-03 

0.4308**1 

17       | 

».11246*,-06 

0.466*370 

18       1 

2.8083*4'-06 

0.46)7718 

14       1 

l.3S69S3*-b6 

0.4421595 

20       1 

6.3561*l'-07 

0.9962122 

21       1 

).1676l5*-37 

0.9911700 

22       1 

1.  3 30**0*  -<,7 

0.4991158 

23       1 

7.39*3***H34 

0.^995728 

2*      1 

).S7259*'-C» 

0.99979)6 

25      1 

1.726136'-U8 

•      8            . 

C      - 

STATE 

2             t             «MU 

P(N<-I> 

1 

Pl»«l  ) 

0.1638197 

1*       1 

0.00021)79035 

0.«)l!S92 

15       1 

0.0031132721 

0.6561621 

16       1 

6.l71*0**-0* 

0.80  50037 

17       | 

3.  )62367*-03 

0.892)99* 

11       1 

l.831919'-95 

0.9*11919 

14       1 

9.9803»»,-3» 

0.9679**1 

20       1 

5-*)7867*-04 

0.982533* 

21       1 

2.)*271»*-36 

9.990*€53 

22       1 

1.61*17*'-G6 

0.99*8162 

23       1 

1.79*5lS*-97 

0.9971737 

2*       1 

*.791*l7*-07 

0.448*612 

25       1 

2.610562*-37 

0.4441616 

26       1 

l.*22Jl)*-37 

0.4445*32 

27       1 

7. 7*91»**-08 

•       8             . 

C      - 

STATE 

2           ■           ""I 

#1  «t<-  1 1 

1 

PIN-I 1 

0.1336713 

16 

0.000279*791 

0.36*1296 

17 

0.00017058*5 

0.378383* 

18 

0.00010*1197 

0.73*7978 

14 

«.39*.0*1*-OS> 

0.8  366240 

20 

).378942*-JS 

0.4000702 

21 

2.)67*30*-05 

0.43*4877 

2i 

l.*»5362*-35 

0.9627610 

2) 

8.8201*8* -J6 

0.4772713 

2* 

5.383507*-O6 

0.986127* 

25 

3.28S902*-9» 

0.4913)27 

26 

2.005598*-06 

0.49*8318 

27 

1.22*l*6'-96 

0.4468*39 

21 

I       7.»7l757*-07 

0.99  807»6 

24 

».<60»96*-07 

0.49822*8 

30 

I       2.78)56S*-07 

0.9442827 

31 

t.698990*-07 

•      8           . 

C     - 

STATE 

2            •            KHO 

fMNX-l I 

I 

*MK«I  ) 

0.10*7044 

14 

1       0.000)539121 

0.2932401 

20 

I       0.0002383763 

0.68  76*89 

21 

I       0.0001621989 

0.6*31621 

22 

1       0.0001192738 

0.7571442 

23 

1        7.*97l63*-3S 

0.93*64*9 

2* 

I       5.947079*-03 

0.8873427 

25 

1        3.*65)19*-04 

0.4235789 

26 

I       2.3S597l*-0* 

0.4*80**7 

27 

1       1.6017*8*  -J* 

0.94*6775 

23 

I        1.038977'-O5 

0.975935* 

24 

|       7.  603601* -06 

9.983*733 

30 

1        3.033*7*'-06 

0.9888444 

31 

1       3.*22096*-9k 

0.942*535 

32 

I       2.32**73*-06 

0.44*849* 

33 

1       1.5«l76l*-96 

0.4969118 

3* 

I       1.073)!8*-46 

0.4976285 

33 

1       7.311221*-37 

0.9983877 

36 

|       *.970663*-O7 

0.4484038 

37 

I        ).37939i'-07 

Pt»<»ll 

1 

0.073783*6 

25       I 

0.22*21*5 

26       1 

0.338381* 

27       1 

0. 532C062 

28       1 

9.6*56931 

29       1 

0.7329*11 

30       1 

0.7987736 

31       1 

0.8*9)17* 

32      1 

0. 8457673 

33       1 

9.9139310 

3*      1 

0.9351507 

35       1 

0.9511389 

36       1 

0.463185) 

37       1 

0.9722617 

31       1 

0.9791003 

34       1 

0.99*2530 

60       1 

0.9891333 

*1       1 

0.991960* 

♦2       1 

0.9432*** 

»)       1 

0.99*9230 

**       1 

3.9961T62 

-■>       1 

0.9971189 

«6        1 

0.4973293 

*7       1 

0.999)6** 

*1       1 

0.993  76/7 

*9       1 

PIN-I  I 

C0003038230 
0.00022I>41»3 
3.000172*778 
0.00012995*6 
4.74l527*-0» 
7.)77**S*-U» 
*.5SS*2*,-OS 
*.1«B18S*-0S 
3.l55618*-05 
2.377623'-05 
1.791*17*-03 
I.  3*  97  70* -Oi 
l.01699»*-3» 
r.66281k'-06 
5.773431 '-*J» 
*. 3 5 00*8* -06 
).27757**-0« 
2.*6«3ll*-06 
1.9a067l'-06 
l.»0143J'-36 
1.036248'-06 
7.451735*-U7 
5.  4  46  530* -0  7 
».5l8166*-07 
l.»042*»*  -07 


P(N<-II 

0.4  44400  2 
0.9999518 
0.944976  7 
0.9999617 
0.99*794*5 
0.444447) 
9.9949987 
0. 99*9991 
0.9949947 
0.9444441 
0.4444499 
0.4499999 
0.4444449 


PINX-I) 

0.4447311 

0.49486** 
0.944  9261 
0.9999597 
0.9999741 
0.9999S3U 
0, 94444)* 
0.494496* 
0.4449410 
0.4444484 
0.499944* 
0.4449496 
C. 4494998 
C.9444499 


Pl-*»ll 

0.4495621 
0.9497327 
0.9998369 
0.499930* 
0.9949)92 
0.999962  4 
0.494477* 
0.9994462 
0.999991  5 
0.99V9'>*3 
0.944446  8 
0.449498U 
0.444449  3 
0.444494) 
0. 9999994 
b.  94494  9  7 


PIN<-I I 

0.99925*0 
0.999*9)3 
0.999653* 
0.9997658 
0.9448*07 
0.4998417 
0.444426* 
0.4444*44 
0.4194660 
0.99997oS 
0.44448*3 
0.  444  7P43 
0.444442  7 
0.4  744451 
0.9449466 
0.9999477 
0.999949* 
0.4449484 
0.4444492 


»|N<-I I 

0.9490  71* 
0.444300* 
0.144*724 
0.4996028 
0.9947907 
0.44977** 
0.4491301 
0.4941720 
0.444433* 
0.44442  7  3 
0.4444*52 
0.4999*87 
0.9999689 
0.9999766 
0.9994823 
0.9999866 
0.9444844 
0.944492* 
0.49949*3 
0.9994937 
0.9  79944  7 
0.94V94/6 
0.9»94992 
0.9949496 
0.4444484 


46 


M      »       i 


K      =      8 


C      =      2 


RHO 

P( DELAY  I 

LIGIVEN  K) 

LQ(  GIVEN  K) 

LQ  FOR  K*l 

RATIO 

0.10 

0.01792109 

0.2013269 

0.301326923 

0.002020202 

0.6568266 

0.20 

0.06526637 

0.4105492 

0.01054924 

0.01666667 

0.6329545 

0.30 

0.1352489 

0.6365210 

0.03652108 

0.05934066 

0.6154479 

0.40 

0.2233869 

0.8917542 

0.09175420 

0.1523809 

0.6021371 

0.50 

0.3265101 

1.197244 

3.1972441 

0.3333333 

0.5917324 

0.55 

0.3829307 

1.380202 

0.2802032 

0.4770609 

0.5873531 

0.60 

0.4422594 

1.593808 

3.3938082 

0.6750000 

0.5834197 

0.65 

0.5042829 

1.851505 

0.5515060 

0.9510822 

0.5798720 

0.70 

0.5688107 

2.175664 

0.7756644 

1.345098 

0.5766602 

0.75 

0.6356713 

2.606503 

1.106503 

1.928571 

0.5737424 

0.80 

0.7047099 

3.224415 

1.624415 

2.844444 

0.5710835 

0.85 

0.7757854 

4.216932 

2.516932 

4.426126 

0.5636536 

0.90 

0.8467694 

6.146583 

4.346582 

7.673684 

0.5664272 

0.95 

0.9235438 

11.82589 

9.925892 

17.58717 

0.564382  3 

0.98 

0.9692222 

28.73332 

26.77332 

47.53494 

0.5632344 

0.99 

0.9845791 

56.86909 

54.88910 

97.51749 

0.5628642 
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COMPUFER    PROGRAM 


PREAMBLE 

PERMANENT    ENTITIES 

THE    SYSTEM   OWNS    A    QUEUE 

EVERY   SERVER    HAS   A    STATE, AN    IOENTI 

TEMPORARY    ENTITIES 

EVERY    JOB    HAS    AN    ARR.TIME    AND    MAY    BELONG   TO    THE    QUEUE 

EVENT    NOTICES    INCLUDE    ARRI VAL,END.OFoSI MULATI ON 

EVERY    DEPARTURE    HAS    A    CODE 

PRIORITY   ORDER    IS    ARRI VAL, DEPARTURE, END. OF. SI MULATI ON 

NORMALLY, MODE    IS    REAL 

TALLY   AV.Q.T    AS    THE    MEAN    AND    VAR.Q   AS    THE    VARIANCE    OF   R 

TALLY    FREQ(0    TO    15    BY    0.25)    AS    THE    HISTOGRAM    OF    R 

TALLY    AVER    AS    THE    MEAN    AND    VARER    AS    THE    VARIANCE    OF    ERLG 

ACCUMULATE    MMM    AS    THE    MEAN    OF    K 

ACCUMULATE    KKK    AS    THE    SUM    OF    MM 

ACCUMULATE    FRES(0    TO    40    BY    1)    AS    THE    HISTOGRAM    OF    K 

DEFINE    FREE    TO   MEAN    0 

DEFINE    BUSY    TO    MEAN    1 

DEFINE    X,Y,XX    AND    YY    AS    REA L, 1-DI MENS IONAL    ARRAYS 

DEFINE    F,Z   AS    REAL , 1-DI MENSIONAL    ARRAYS 

DEFINE    NN,KK,STP,EPSTP,SToN,EP    AS    INTEGER    VARIABLES 

DEFINE    U,ERLG,T0T1,AVG,MM    AS    REAL    VARIABLES 

DEFINE    STS1,STS2,STS3    AS    INTEGER    VARIA6LES 

DEFINE    SD1,SD2,SD3    AS    INTEGER   VARIABLES 

DEFINE    FLAG    AS    AN    INTEGER    VARIABLE 

DEFINE    LAMBDA,MU,S,T0ToQoTIME,R,SUMSQR,D,BETA,ST,STT,STARToT 

AS    REAL    VARIABLES 
DEFINE    CUSTOMER, J, K, DELAYER, AR,M,N, ALPHAS A2E.N0, I, EFLAG 

AS    INTEGER   VARIABLES 
END 


MAIN 

READ   CASE. NO 

READ    J 

• •    J    IS    THE    NUMBER    OF    SERVERS 

READ    ALPHA 

»■  ALPHA  IS  SHAPE  PARAMETER  OF  THE  SERVICE  TIME  DISTRIBUTION 

READ  BETA 

••  BETA  IS  SCALE  PARAMETER  OF  THE  SERVICE  TIME  DISTRIBUTION 

READ  LAMBDA 

■  ■  LAMBDA  IS  ARRIVAL  RATE 

READ  EPSTP 

'•  EPSTP  IS  TOTAL  NUMBER  OF  EPOCHS  FOR  THIS  SIMULATION  RUN 

READ  FLAG 

'•  FLAG  INDICATES  WHICH  PAIR  OF  ANTITHETIC  VARIATES  WILL  BE 

■■  USED 

READ  SC1,SD2,SD3 

READ    NN.KK 

•  •    NN    IS    DIMENSION    OF    ARRAY    F 

11     KK    IS    NUMBER    OF    ELEMENTS    NE    1    IN    ARRAY    F 

RESERVE    F(*)    AS    NN    AND    Z(*)    AS    NN 

••     READ    CDF    TABLE    OF    THE    CHI-SQR    DISTN.    WITH    2*ALPHA 

• •  DEGREES    OF    FREEDOM 

FOR    1=1    TO    KK,DO 

READ    F(I) 

LET    F(I)=1.0-F(I) 
LOOP 
FOR    1=1    TO    KK,DO 

READ    Z(I) 
LOOP 
FOR    I=KK+1    TO    NN,DO 

LET    F(I)  =  1.0 
LuQD 
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LET    N.SERVER=J 

CREATE    EACH   SERVER 

LET    K=0 

LET    STP=0 

LET    STS1=SEED.V(SD1) 

LET    STS2=SEED«V(SD2) 

LET    STS3=SEED.V(SD3) 

START    NEW    PAGE 

SCHEDULE    AN    ARRIVAL    NOW 

START    SIMULATION 

STOP 

END 


EVENT  ARRIVAL  SAVING  THE  EVENT  NOTICE 

DEFINE  JJ  AS  AN  INTEGER  VARIABLE 

ADD  1  TO  AR 

LET  T=TIME.V*1440 

IF  STP  EQ  1 

DESTROY  THIS  ARRIVAL 

RETURN 

OTHERWISE 

ADD  1  TO  K 

••  GENERATE  INTERARRIVAL  TIME 

LET  U=UNIFORM.F(OoO,1.0fSDl) 

IF  FLAG  EQ  1 

LET  L=loO-U 

ALWAYS 

LET  AR.TIME=-LOG.E.F(U)/LAMBDA 

RESCHEDULE  THIS  ARRIVAL  IN  AR.TIME  MINUTES 

LET  JJ=0 

FOR    EVERY    SERVER, DO 

IF    STATE(SERVER)=BUSY 

ADD    1   TO    JJ 

ALWAYS 

LOOP 

IF    JJ=0 

LET    EFLAG^l 

CALL    EPOCH 

GO   TO    AA 
ALWAYS 
LET    J1=J-1 
IF    JJ    EQ    Jl 

LET    MM=1.0 
ALWAYS 

•  AA  ■ 

•NEXT*       FOR    EVERY   SERVER, DO 
IF    STATE(SERVER)=FREE 

LET    STATE(SERVER)=BUSY 
LET    R=OoO 
•ADD1       ADD    1    TO    CUSTOMER 

•  ■    GENERATE    SERVICE    TIME 
LET    U=UNIF0RM.F(0o0,l«0,SD2) 
IF    FLAG    EQ    i 

LET    U=1.0-U 

ALWAYS 

CALL    ERLNG 

LET    DEP.TIME=ERLG 

LET    IDENTI(SSRVER)=DEP„TIME 

SCHEDULE    A    DEPARTURE    GIVEN    DEP.TIME    IN    DEP.TIME    MINUTES 

RETURN 

OTHERWISE 

LOOP 

CREATE    A    JOB 

LET    ARRoTIME( JOB )=TIME» V*1440 

FILE    THIS    JOB    IN   THE    QUEUE 

RETURN 

END 
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EVENT  DEPARTURE  SAVING  THE  EVENT  NOTICE 

DEFINE  JJ  AS  AN  INTEGER  VARIABLE 

LET  T=TIME.V*1440 

IF  STP  EQ  1 

DESTROY  THIS  DEPARTURE 

RETURN 

OTHERWISE 

SUBTRACT  1  FROM  K 

LET  JJ=0 

FOR    EACH   SERVER, DO 

IF    STATE(SERVER)=BUSY 
ADD    1    TO    JJ 

ALWAYS 
LOOP 
IF    JJ    EC    1 

LET    EFLAG=0 

CALL    EPOCH 

GO   TO    AA 
OTHERWISE 
IF    JJ    EQ    J 

IF    N.CUEUE    EQ    0 
LET    MM^OoO 

ALWAYS 
ALWAYS 
•  AA ' 
FOR    EVERY   SERVER, DO 

I F    I  DENT  I ( SERVER ) =CODE (DEPARTURE  I 
IF    THE    QUEUE    IS    EMPTY 

LET    STATE(SERVER)=FREE 
DESTROY    THE    DEPARTURE 
RETURN 
OTHERWISE 

LET    AoT=ARR.TIME(FoQUEUE) 
REMOVE    FIRST    JOB    FROM    THE    QUEUE 
DESTROY    THIS    JOB 
LET    R=T-AoT 
ADD    1    TO    CUSTOMER 
ADD    I   TC    DELAYER 
ADD    R    TO   TOToQ«TIME 
ADD    R**2   TO    SUMSQR 
• •     GENERATE    SERVICE    TIME 
LET    U=UNIFORM.F(0.0,1.0,SD3) 
IF    FLAG    EQ    1 
LET    L=1«0-U 
ALWAYS 
CALL    ERLNG 
LET    DEPoTIME=ERLG 
LET    IDENTI(SERVER)=DEPoTIM5 

SCHEDULE    A    OEPARTURE    GIVEN    DEP.TIME    IN    OEP.TIME    MINUTES 
RETURN 
OTHERWISE 
LOOP 
END 


EVENT    END.OFoSIMULATION 

DEFINE    FR, LUCKY    AS    INTEGER    VARIABLES 

DEFINE    DELI    AS    AN    INTEGER    VARIABLE 

RESERVE    X(*)    AS    60, Y(*)    AS    60    AND    YY(*)    AS    60 

RESERVE    XX(*)    AS    60 

LET    S=TIMEaV*1440 

START    NEW    PAGE 

LET    MU=1.0/(ALPHA*BETA) 

•  •  MU  IS  SERVICE  RATE 
LET  G=LAMBDA/MU 

•  •  G  IS  OFFERED  LOAD 
LET  RO=G/J 
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••    RG    IS   TRAFFIC    INTENSITY 

»■     COMPUTE    AVG. WAIT. TIME    FOR    EQUIVALENT    M/M/C    QUEUE 

LET    JFAC=1 

FOR    1=1    TO    J-l.DO 

LET    JFAC=JFAC*(I+1> 
LOOP 

LET  C=((J*RO)**J)/(JFAC*(loO-RO) ) 
LET  B=1.0 
FOR  K=l  TO  J-i,DO 
LET  KFAC=1 

FOR  L=l  TO  K-1,D0 

LET  KFAC=KFAC*(L+1) 

LOOP 
ADD    <(J*RO)**K)/KFAC   TO    B 
LOOP 

LET    PCNE=B+C 
LET    PNGT=1.0/(B+C) 

LET    E=PNOT*{ (J*RO)**J)/(J*JFAC*(( 1<»0-R0)**2  )*MU) 
LET    VSQR=loO/ ALPHA 
LET    V=SQRT.F( VSQR) 
LET    W=(E/2oO)*(i«Q+VSQR) 
LET    PMMC=E*J*MU*(l«0-RO) 
PRINT    1    LINE    WITH   CASE. NO    AND    FLAG    AS     FOLLOWS 

CASE   NO=****  FLAG=* 

SKIP    2    LINES 

PRINT    2    LINES    AS    FOLLOWS 
INPUT    INFORMATION: 

5RTPTCTRK 

PRINT  5  LINES  WITH  STS1 ,SEEDo V( SD1 ) ,STS2,SEED. V (SD2 ) i 
S7S3,ScED.V(SD3)  AS  FOLLOWS 

STARTING  SEEDS:  FINAL  SEEDS: 

I'~i~*"*5*55***"*"  #*"*"*"*"***"** 

I   a     ^*  " P  ^P  ^H  ^P  "^  T*  ^T  ^H  ^J*  ^^  *?  ^P  ^F  ^P  ^P  ^P  ^P  *f*   ^^ 

III.   **********  ********** 

SKIP  2  LINES 

PRINT  1  LINE  WITH  ALPHA  AMD  BETA  AS  FOLLOWS 

ALPHA=***      BETA=**.**** 
SKIP  2  LINES 
PRINT  1  LINE  WITH  LAMBDA  AS  FOLLOWS 

LAMBDA=***„**** 
SKIP  2  LINES 
PRINT  1  LINE  WITH  MU  AS  FOLLOWS 

MU= ***„**** 
SKIP  2  LINES 

PRINT    1    LINE    WITH    RO    AS    FOLLOWS 
TRAFFIC    INTENSITY, RO=*.**** 
SKIP    2    LINES 

PRINT    1    LINE    WITH    J    AS    FOLLOWS 
NUMBER    OF   SERVERS=** 
SKIP    2    LINES 

PRINT    1    LINE    WITH    G    AS    FOLLOWS 
OFFERED    LOAD=**.**** 
SKIP    2    LINES 

PRINT    1    LINE    WITH    V    AS    FOLLOWS 
COEFFICIENT   OF    VARI ATION=**.**** 
START   NEW    PAGE 
IF    EFLAG=1 

ADD    S-START.T    TO   SUMP 
ALWAYS 

LET    LUCKY=CUSTGM5R-DELAYER 
LET    PD=DELAYER/CUSTOMER 
PRINT    2    LINES    AS    FOLLOWS 
OUTPUT    INFORMATION: 

SKlp-rTTFiEs^ 

PRINT  1  LINE  WITH  AV.Q.T  AS  FOLLOWS 

MEAN  QUEUEING  TIME  FROM  SIMULAT ION=***o**** 

SKIP  2  LINES 

PRINT  1  LINE  WITH  W  AS  FOLLOWS 

COMPUTED  MEAN  QUEUEING  TIME =***.****  MINUTES 
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SKIP    2    LINES 

PRINT    1    LINE    WITH   W-AVoQoT    AS    FOLLOWS 

DIFFERENCE****.**** 

SKIP  2  LINES 

PRINT  1  LINE  WITH  VAR.Q  AS  FOLLOWS 

VARIANCE  OF  QUEUEING-TI ME=***.**** 

SKIP  2  LINES 

PRINT  1  LINE  WITH  MMM  AS  FOLLOWS 

AVERAGE  QUEUE  LENGTH***.**** 

SKIP  2  LINES 

PRINT  1  LINE  WITH  CUSTOMER  AS  FOLLOWS 

TOTAL  NUMBER  SERVED  IN  SIMULATION  PERIOD****** 

SKIP  2  LINES 

PRINT  1  LINE  WITH  DELAYER  AS  FOLLOWS 

TOTAL  NUMBER  OF  CUSTOMERS  WITH  DELAY  GREATER  THAN  ZERO=***** 

LET  PERIOD=S-STT 

SKIP  2  LINES  I 

PRINT  1  LINE  WITH  PD  AS  FOLLOWS 

PROB(W>0)=*. ****** 

SKIP  2  LINES 

PRINT  1  LINE  WITH  PMMC  AS  FOLLOWS 

PMMC=*. ****** 

SKIP  2  LINES 

PRINT  1  LINE  WITH  PERIOD  AS  FOLLOWS 

SIMULATION  TIME*******.****   MINUTES 

SKIP  2  LINES 

PRINT  1  LINE  WITH  KKK*1440  AS  FOLLOWS 

TOTAL  BUSY  PERIOD*********.****   MINUTES 

FOR  1*1  TO  60, DO 

LET  XX(I)*FREQ(I) 
LOOP 

SUBTRACT  LUCKY  FROM  XX( 1 ) 
»■  COMPOTE  AND  PRINT  F  SUB  I 
FOR  1=1  TO  N,DO 
ADD  XX(I)  TO  FR 
LET  Yd  )  =  1.0-(FR/DELAYER) 
IF  FLAG  EQ  0 
IF  FR  EC  DELAYER 
LET  N*I-1 
GO  TO  REG 
OTHERWISE 
ALWAYS 

LET  Y(I)=LOG.E.F(Y(I)) 
LET  X(I  )=IM.O 
LOOP 
•REG* 

LET  CUMFREQ=0 
START  NEW  PAGE 

■■  COMPUTE  AND  PRINT  STATE  PROBABILITIES 
PRINT  2  LINES  AS  FOLLOWS 
I  X(I)  FREQ(I)  CUMFREQ 


PUR  1=1  TO  N,D3 

ADD  XX(I)  TO  CUMFREQ 

SKIP  1  LINE 
PRINT  1  LINE  WITH  I,X(I),XX(I)  AND  CUMFREQ  AS  FOLLOWS 
**  **.***  *********  ********* 

LOOP 

IF  FLAG  EQ  1 

CALL  REGRESSION 
ALWAYS 

START  NEW  PAGE 
PRINT  2  LINES  AS  FOLLOWS 
N  T(N)  P(N) 

FOR"  1=1  TO  407D"0 

PRINT  1  LINE  WITH  I , FRES ( I ) *1440.0  AND  ( FRES ( I )*1440. ) /PERIO 

AS  FOLLOWS 
**  *****,****  *.******** 

LOOP 

RETURN 

END 
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ROUTINE  EPOCH 

DEFINE  FI.N  AS  AN  INTEGER  VARIABLE 

LET  T=TIME.V*1440 

IF  EFLAG  EQUALS  1 

•«  EPOCH  STARTS 

LET  START. T=T 

LET    ST.N=CUSTOMER 

LET    TOTl=TOToQ.TIME 

ADD    1    TO    EP 

RETURN 

OTHERWISE  v 

• ■     EPOCH    ENDS 

LET    FINISH.T  =  T 

LET    FI.N=CUSTOMER 

LET  TOT2=TOT.Q.TIME 

LET  AVG=(TOT2-TOTl)/< ( F I.N-ST.N )*1.0 ) 

LET    PERIOD=FINISHoT-START«T 

LET    CUM.AVG=T0T2/FI.N 

SKID    1    LINE 

LIST    EP,START.T,FINISHoT,ST.N,FI.N,T0Tl,TQT2 

LIST    PERIOD,TOT2-TOTltFIoN-ST«N,AVG,CUM.AVG 

•«    ChECK    FOR    END   OF    SIMULATION 

IF    EP    EQ    EPSTP 

LET    STP=1 

SCHEDULE    AN    END. OF. SI MULATION    NOW 
ALWAYS 
RETURN 
END 


ROUTINE    ERLNG 

DEFINE    PL,LH,RH    AS    INTEGER    VARIABLES 

IF    U    LT    F(l) 

LET    CHSQR=ZU)*(U/FU)  ) 
GO    TO    ERLANG 
ALWAYS 

••     START    BISECTION    SEARCH 
IF    U    GT    F(KK) 
LET    CHSQR=Z<KK)+(Z<KK)-Z(KK-l))*(U-F(KK)) 

GO    TO    ERLANG 
ALWAYS 
LET    LH=0 
LET    RH=NN 
•AA'    LET    RL=(LH+RH)/2 

IF  U  LT  F(RL) 
LET  RH=RL 
GO  TO  BB 
ALWAYS 

LET  Lh=RL 
•BB'  IF  RH=LH+1 

GO  TO  QUIT 
ALWAYS 
GO  TO  AA 
•QUIT* 

11  DO  LINEAR  INTERPOLATION 

LET  CHSQR=Z(LHH-((U-F(LH))*<Z(RH)-Z(LH)  )  )/ (  F(  RH) -F(  LH)  I 
•ERLANG1 

LET  EPLG=(BETA*CHSQR)/2.0 
RETURN 
END 
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